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SECTION  1 
INTRODUCTION 


This  report  documents  numerical  algorithms,  numerical  examples,  and  parametric 
computations  for  crack  propagation  in  two  and  three  dimensions.  Volume  1  of  this  report 
presents  background  information,  experimental  data,  analysis  of  the  data,  and  computational 
approaches  and  results  [1].  The  numerical  algorithms  described  in  this  report  have  been 
incorporated  into  the  1999  version  of  the  EPIC  code  [2]. 

The  objective  of  this  work  was  to  develop  the  capability  to  perform  3D  crack  propagation 
computations  for  penetrator  cases  during  impact  and  penetration  into  hard  concrete  targets.  The 
first  step  was  to  develop  a  new  2D  algorithm  that  did  not  require  rezoning  (as  rezoning  would  be 
very  complex  in  3D  geometry),  the  second  step  was  to  extend  the  new  technique  to  3D  geometry, 
and  the  third  step  was  to  apply  the  3D  algorithm  to  actual  hard  target  penetrator  impact  events. 


Some  background  on  computational  approaches  for  crack  propagation  is  as  follows: 

Finite  element  modeling  of  fracture  propagation  has  traditionally  been  achieved  by 
splitting  nodes  along  a  crack  path  and  allowing  the  adjacent  elements  to  separate.  If  there  are  no 
symmetry  assumptions  and  the  crack  is  allowed  to  propagate  in  any  direction,  remeshing 
algorithms  become  necessary;  examples  of  such  algorithms  appear  in  [3,  4,  5,  6].  If  a  plane  of 
symmetry  is  assumed  —  as  is  often  the  case  when  the  study  does  not  focus  on  the  direction  of 
crack  propagation  —  the  node-splitting  technique  reduces  to  a  “node-release”  technique  in  which 
the  propagating  crack  is  modeled  by  the  release  of  restrained  nodes  along  the  line  of  symmetry. 
Examples  of  this  type  appear  in  [7,  8,  9].  In  the  node-release  technique,  external  forces  are 
typically  applied  to  the  nodes  immediately  after  they  are  released.  These  forces  are  then  reduced 
to  zero  as  the  crack  tip  traverses  the  element  edge,  reducing  the  waves  generated  by  a  sudden 
release  of  restraint  and  extracting  energy.  For  both  the  node-splitting  and  node-release 
techniques,  a  fracture  criterion  (in  conjunction  with  an  energy  integral  or  the  crack-tip  opening 
angle)  must  be  employed  to  determine  the  speed  of  propagation.  For  the  node-splitting  technique 
(in  general  directions),  a  criterion  to  determine  the  direction  of  propagation  must  also  be 
employed. 
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As  an  alternative  to  the  node-splitting  technique,  micromechanically  based  approaches 
have  been  proposed  in  which  the  speed  and  direction  of  propagation  are  determined  directly  from 
the  constitutive  models  via  material  failure.  Examples  of  these  approaches  appear  in  [10,  11]. 
To  elaborate  on  a  specific  example,  Needleman  and  Tvergaard  [10]  use  an  elastic-viscoplastic 
constitutive  relation  for  a  porous  plastic  solid,  which  models  the  nucleation  and  growth  of  voids 
to  coalescence.  They  demonstrate  a  loss  of  convergence  of  a  dynamic  J-integral  with  mesh 
refinement  when  a  void-size  length  scale  is  omitted  from  the  constitutive  relation.  An  advantage 
of  the  micromechanically  based  approaches  is  the  relative  ease  with  which  the  crack  model  is 
propagated  through  the  mesh;  remeshing  and  modifying  contact  surfaces  are  not  necessary. 

Dynamic  fracture  has  also  been  modeled  by  means  other  than  finite  elements,  e.g.,  the 
boundary  element  (BE)  method  and  the  element-free  Galerkin  (EFG)  method.  Belytschko  et  al 
[12,  13,  14]  use  the  EFG  method,  in  which  the  spatial  discretization  consists  only  of  nodes,  and 
shape  functions  are  generated  by  a  moving  least-squares  approximation  of  the  nodal  data.  This 
method  allows  the  advancing  crack  to  be  modeled  by  simply  modifying  the  nodal  connectivities. 
Energy  integrals  are  used  as  fracture  parameters. 

In  this  report,  a  finite  element  algorithm  is  presented  for  dynamic  crack  calculations  in 
general  directions.  The  algorithm  avoids  node  splitting  and  the  associated  remeshing  and 
redefinition  of  contact  surfaces,  and  does  not  rely  on  failure  mechanisms  embedded  in  the 
constitutive  model.  The  T*  energy  integral  is  employed  as  the  dynamic  fracture  parameter, 
guiding  the  crack  tip  through  the  mesh.  Elements  through  which  the  crack  tip  passes  lose  the 
ability  to  sustain  deviatoric  and  tensile  volumetric  stresses,  and  their  interfaces  with  the  rest  of 
the  mesh  model  fractured  surfaces  that  can  sustain  only  compressive  normal  tractions  (when  the 
crack  is  closed).  This  technique  is  therefore  referred  to  as  the  element-failure  algorithm. 
Discussions  of  crack-tip  energy  integrals  leading  to  the  development  of  T*  appear  in  Atluri  et  al 
[15]  and  Nakamura  et  al  [16].  The  development  of  T*  appears  in  Atluri  et  al  [17, 18]. 

The  following  three  sections  provide  a  description  of  the  2D  algorithm,  a  description  of 
the  3D  algorithm,  and  a  computational  parametric  study  for  a  hard  target  penetrator. 
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SECTION  2 

2D  ALGORITHMS 


2.1  THE  2D  ELEMENT-FAILURE  ALGORITHM 

Figure  1  demonstrates  a  typical  remeshing  technique  used  with  a  node  splitting  algorithm. 
The  region  of  the  mesh  near  the  crack  tip  is  remeshed  as  the  tip  moves  through  the  mesh, 
allowing  propagation  in  any  direction.  Remeshing  entails  optimizing  the  new  mesh,  modifying 
the  nodal  connectivites  accordingly,  and  mapping  state  variables  from  the  old  mesh  to  the  new 
one.  The  new  crack  surfaces  are  then  incorporated  into  the  contact  algorithm. 


In  the  element-failure  algorithm  of  Figure  2,  a  criterion  for  the  direction  of  propagation  is 
invoked  each  time  the  crack  tip  enters  a  new  element.  (The  tip  is  restricted  to  traveling  in  a 
straight  line  through  the  element.)  After  the  direction  has  been  determined,  the  fracture 
parameter  T*,  when  combined  with  a  fracture  criterion,  determines  the  speed  of  propagation.  T 
is  recalculated  at  each  timestep,  and  the  tip  speed  is  subsequently  updated.  The  elements  through 
which  the  tip  passes  lose  the  ability  to  sustain  deviatoric  and  tensile  volumetric  stresses,  i.e.,  they 
are  failed.  This  technique  requires  no  remeshing,  and  none  of  the  associated  optimization  and 
mapping.  The  failed  elements  (shaded  in  Figure  2)  can  sustain  only  positive  pressure,  and  this 
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Failed  elements  (shaded)  Cracktip  path 


Figure  2.  Crack  Path  Model  of  the  Element-failure  Algorithm  [02.jr103794.tifi 

typically  occurs  under  negative  volumetric  strains.  An  edge  between  a  failed  and  an  unfailed 
element  therefore  acts  as  a  traction-free  crack  face  when  the  volumetric  strain  of  the  failed 
element  is  positive  (opened  crack),  and  a  pressure-transmitting  crack  face  when  the  volumetric 
strain  of  the  failed  element  is  negative  (closed  crack).  Crack-face  surfaces  are  not  incorporated 
into  the  contact  algorithm. 

If  a  crack  path  falls  along  a  plane  of  symmetry,  a  node-splitting  algorithm  may  be  reduced 
to  a  node-release  algorithm.  The  left  side  of  Figure  3  shows  how,  in  a  node-release  algorithm, 
the  advancing  crack  is  modeled  by  the  release  of  nodes  which  were  restrained  in  the  x-direction 
along  the  line  of  symmetry.  The  sudden  release  of  these  nodes  generates  release  waves  (noise), 
an  effect  that  is  often  avoided  [7,  8,  9]  by  replacing  the  x-restraint  with  an  externally  applied 
nodal  force.  This  force,  denoted  by  P  on  the  left  side  of  Figure  3,  is  effectively  equivalent  to  the 
restraint,  but  it  is  subsequently  reduced  to  zero  as  the  crack  tip  moves  toward  the  next  node  along 
the  line  of  symmetry.  This  technique  is  used  for  the  node-release  calculations  reported  in 
subsection  2.3.1. 
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Current  tip  element 

iz 


►  x 


Line  of 

symmetry  Cracktip  path 

Figure  3.  Release  Forces  in  Node-release  Algorithm  (left,  P)  and  Element-failure 
Algorithm  (right,  Pxi,  Pzi, ....  PZ3>  [03_n03794.tif] 

For  the  element-failure  algorithm,  the  right  side  of  Figure  3  shows  the  manner  in  which 
release  forces  are  used  to  avoid  the  release  waves  associated  with  the  sudden  failure  of  the 
current  crack-tip  element.  When  the  crack  tip  first  enters  the  element,  the  element  is  failed  and 
external  forces  (Pxi,  Pzi,  •••,  Pz3  in  Figure  3)  are  applied  to  all  the  nodes  of  the  failed  element. 
The  external  forces  are  equivalent  to  the  internal  forces  on  those  nodes  due  to  the  stresses  in  the 
element  before  it  was  failed.  As  in  the  node-release  algorithm,  the  release  forces  are 
subsequently  reduced  to  zero  as  the  crack  tip  moves  through  the  element.  They  vary  linearly  to 
zero  as  the  tip  reaches  a  distance  of  7  2  Aeie  from  its  location  when  the  forces  were  initially 

applied.  (Aeie  =  current  area  of  the  element.)  Shorter  values  of  this  distance  were  found  to 
generate  significantly  more  noise  in  the  solution. 

2.2  THE  T*  FRACTURE  PARAMETER 

As  a  measure  of  the  energy  release  associated  with  a  unit  area  of  crack  surface  extension 
in  dynamic  fracture,  T*  [17, 18]  is  expressed  as, 

T*=«sU(w+Tk-w>*  (1) 

where  T*  is  a  line  contour  around  the  crack  tip.  W  and  T  represent  the  internal  energy  density 
and  kinetic  energy  density,  respectively;  ay,  the  ij-th  component  of  the  stress  tensor;  nj,  the  j-th 
component  of  the  unit  normal  vector  to  T*;  np,  the  component  of  the  unit  normal  vector  in  the 
direction  of  crack  propagation;  ub  the  i-th  component  of  the  displacement  vector;  uijP,  the 
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derivative  of  the  i-th  component  of  the  displacement  vector  in  the  direction  of  crack  propagation. 
(Note  p  is  not  a  free  subscript.) 

The  largest  errors  in  a  numerical  solution  occur  in  the  region  of  the  crack  tip,  where  field 
quantities  vary  most  rapidly.  As  a  result,  it  is  desirable  to  equate  the  limiting  integral  of 
Equation  1  to  area  and  contour  integrals  by  the  divergence  theorem  [7], 

T*=ri(w+Tk-<wJdr  -JJw+rirP^.p-^kjU'*0  <2> 

where  T  represents  a  large  crack-tip  contour  and  is  the  area  between  T  and  T*.  A  subscript 
following  a  comma  represents  the  direction  of  a  spatial  derivative. 

For  the  computations  of  this  study,  meshes  were  constructed  of  constant-strain  triangles 
in  the  “crossed”  configuration.  Figure  4  shows  how  T*,  T  and  Q  can  be  arranged  around  a  crack 
tip  and  its  associated  failed  elements  in  a  crossed-triangle  mesh.  While  the  definition  of  T*  in 
Equation  1  is  the  near-field  integral  around  T*,  its  calculation  is  achieved  by  the  far-field  integral 
around  T  and  the  area  integral  over  Q.  The  segments  of  T*  in  Figure  4  are  not  used  in  this 
calculation. 
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Figure  4.  Domains  of  Integration  in  the  Calculation  of  T*  [04_Tio3794tm 
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For  each  triangular  element,  the  displacement  is  interpolated  linearly,  resulting  in 
constant  displacement  derivatives,  stresses,  densities,  and  internal  energy  densities.  Numerical 
integration  of  the  integral  over  T  in  Equation  2  is  achieved  by  the  trapezoidal  rule  applied  to  each 
of  the  boundary  segments  in  counterclockwise  order.  This  scheme  requires  nodal  values  of  the 
integrand,  which  includes  the  internal  and  kinetic  energy  densities,  the  stresses,  and  the 
displacement  derivatives.  To  generate  nodal  values  of  these  quantities,  values  from  adjoining 
elements  are  averaged  at  the  nodes. 

The  integral  over  £2  in  Equation  2  is  integrated  by  employing  one-point  quadrature  over 
each  triangular  element  in  Q.  For  constant-strain  triangles,  the  second-order  displacement 
derivatives  must  be  smoothed  before  application  of  the  quadrature  scheme  (because  they  are  zero 
in  the  interiors  of  the  elements  and  infinite  at  the  element  boundaries).  Smoothing  is  achieved  by 
averaging  values  of  the  first  displacement  derivatives  from  the  adjoining  elements  at  each  node, 
and  then  interpolating  those  nodal  values  linearly  over  each  element  to  render  constant  second 
displacement  derivatives  in  each  element.  A  similar  approach  is  used  for  the  energy-density 
derivatives.  Dexter  and  O’Donoghue  [8]  also  use  a  node-averaging  technique  to  smooth  the 
internal  energy  and  strain  derivative  fields,  and  they  report  greatly  improved  accuracy  in  the 
computed  values  of  T*. 


Since  the  limit  in  Equation  1  cannot  be  taken  in  a  numerical  solution,  it  is  necessary  to 
use  an  inner  contour,  T*,  of  finite  size.  The  use  of  the  divergence  theorem  to  generate  Equation 
2  is  motivated  by  the  greater  accuracy  of  the  finite  element  solution  in  the  far-field  region  than  in 
the  near-field  region. 

It  has  been  reported  [8]  that  a  lack  of  convergence  with  mesh  refinement  occurs  for 
solutions  of  ductile  fracture  if  T*  is  taken  at  the  limit  value  of  zero.  In  subsection  2.3.2,  it  will  be 
shown  that  a  finite  value  of  T*  is  necessary  for  convergence.  It  is  also  recommended  [7,  8]  that 
all  elements  that  have  been  inside  T*  be  excluded  from  Q  in  subsequent  evaluations  of  T*  due  to 
the  lack  of  accuracy  of  the  numerical  solution  in  the  crack  tip  region  and  the  continuation  of 
those  errors  when  path-dependent  materials  (plasticity)  are  unloaded.  The  result  is  a  zone  of 
elements  that  are  not  included  in  Q  forming  a  wake  behind  the  crack  tip.  This  exclusion  zone  is 
depicted  in  Figure  4,  and  its  effect  is  a  modification  of  the  shape  of  T*.  In  this  report,  the  size  of 
the  exclusion  zone  will  be  reported  by  a  radius,  and  all  elements  whose  centroids  have  been,  at 
any  time,  within  that  radius  of  the  crack  tip  will  be  excluded.  In  a  similar  manner,  the  size  of  the 
outer  contour  T  will  be  reported  by  a  radius,  and  the  edges  of  the  elements  whose  centroids 
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currently  lie  within  that  radius  of  the  crack  tip  and  outside  the  exclusion  zone  will  determine  D 
The  domain  Q  is  defined  by  the  region  contained  inside  the  closed  contour  T  +  T*. 

2.3  NUMERICAL  RESULTS 

Both  the  2D  and  3D  algorithms  have  been  incorporated  into  the  1999  version  of  the  EPIC 
Code  [2],  which  is  based  on  an  explicit  dynamic  displacement-based  finite  element  formulation 
[19].  The  equations  of  motion  are  integrated  by  central  differences  and  the  coupled  equations  of 
crack  propagation  are  integrated  by  the  forward-Euler  scheme.  All  2D  calculations  were 
performed  in  plane  strain. 

For  the  element-failure  calculations,  no  assumptions  are  made  on  the  direction  of 
propagation  due  to  symmetry.  The  direction  may  be  chosen  perpendicular  to  the  direction  of  the 
maximum  normal  stress  in  the  region  of  the  crack  tip,  implying  mode  one  fracture,  or  as  that 
which  produces  the  largest  value  of  T*,  allowing  mixed-mode  fracture.  Because  this  formulation 
does  not  require  the  generation  of  a  stiffness  matrix,  the  nonlinearities  are  treated  in  an  explicit 
manner,  and  many  of  the  difficulties  associated  with  zeroing  stiffness  terms  are  avoided. 

2.3.1  Comparison  of  Element  Failure  and  Node  Splitting 

The  primary  advantage  of  the  element-failure  method  is  the  simplicity  of  implementation 
—  there  is  no  need  for  remeshing  or  redefining  contact  surfaces.  The  accuracy  of  the  method, 
however,  can  be  gauged  by  a  comparison  of  results  to  the  prevalent  node-splitting  method.  The 
sample  problem  used  for  comparison  is  the  coupled  pressure-bar  specimen  of  Figure  5  under 
dynamic  loading.  The  material  model  is  elastic/perfectly  plastic  with  the  von  Mises  flow  rule 
a/^T  =  ay  ■  The  material  parameters  are  intended  to  represent  a  steel:  elastic  modulus,  E  =  200 

GPa;  Poisson’s  ratio,  v  =  0.3;  and  uniaxial  yield  stress,  ay  =  1.1  GPa.  The  fracture  response  is 
prescribed  by  the  resistance,  R(vc),  a  function  of  tip  velocity,  and  it  is  plotted  on  the  right  side  of 
Figure  5.  This  resistance  exhibits  a  very  strong  dependence  on  crack  speed,  and  it  was  selected 
to  produce  substantial  crack-tip  plasticity.  For  values  of  T*  less  than  the  threshold  value  of 
R(vc  =  0)  =  20  kN/m,  the  crack  tip  velocity  vc  is  zero.  For  values  of  T*  above  the  threshold,  the 
resistance  increases  linearly  and  T*  =  R(vc)  determines  the  tip  velocity.  The  loading  is  a 
prescribed  velocity  that  increases  linearly  with  time, 

v*  =  <M  (3) 
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where  av  =  1 10,000  m/s2;  see  Figure  5. 


Figure  5.  Coupled  Pressure-bar  Specimen  and  Fracture  Response  Function 

[05_T103794.tif] 

Symmetry  along  the  crack  path  allows  for  use  of  the  node-release  algorithm.  Its  mesh  is 
shown  on  the  left  side  of  Figure  6,  and  contains  1419  nodes  and  2720  elements.  For  the  element- 
failure  algorithm,  a  full  mesh  of  2862  nodes  and  5536  elements  is  used,  as  shown  on  the  right 
side  of  Figure  6.  For  both  algorithms,  the  radius  off  is  1.5  cm  and  the  radius  of  T*  is  0.5  cm. 


Figure  6.  Meshes  for  Comparison  of  Node-release  (left,  2720  elements)  and  Element- 
failure  (right,  5536  elements)  Methods  [06a/b_T1 03794.tif] 
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Figure  7  contains  plots  of  T*  and  the  z-coordinate  of  the  crack  tip  as  functions  of  time  for 
both  the  node-splitting  and  element-failure  methods.  Close  agreement  indicates  similar  accuracy 
of  the  two  methods.  Figure  8  is  a  plot  of  the  mesh  at  250  ps,  with  the  failed  elements  shaded. 
This  plot  verifies  adequate  performance  of  the  propagation-direction  criterion. 


Figure  7.  Comparison  of  Results  for  Node-release  and  Element-failure  Methods 

[07a/b_T 1 03794.tif] 
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Figure  8.  Element-failure  Mesh  at  250  ps  (failed  elements  shaded)  [08.Tio3794.tin 


The  large  drop  in  the  value  of  T*  in  Figure  7  around  210  ps  is  a  result  of  the  crack 
approaching  the  top  surface  of  the  specimen.  When  the  distance  from  the  crack  tip  to  a  surface  is 
shorter  than  the  radius  of  T*,  as  depicted  in  Figure  9,  a  segment  of  T*  will  coincide  with  a 
segment  of  T.  The  effect  of  this  coincident  segment  is  the  loss  of  the  contribution  of  that 
segment  to  the  computed  value  of  T*,  and  a  corresponding  reduction  in  the  speed  of  the  crack  tip, 
which  can  be  seen  as  a  change  in  slope  of  the  plot  of  the  z-coordinate  of  the  crack  tip  versus  time 
in  Figure  7. 


Figure  9.  Domains  of  Integration  Near  a  Surface  [09.Ti03794.tiq 
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2.3.2  Convergence  of  T* 


In  this  section,  the  sample  problem  of  subsection  2.3.2  is  used  to  study  the  convergence 
characteristics  of  the  T*  fracture  parameter  with  mesh  refinement.  It  has  been  shown  [7,  8]  that 
the  node-release  algorithm  converges  for  a  constant  radius  of  T*.  Brust  et  al  [7]  conclude  that 
the  radius  of  T*  should  be  chosen  small  enough  so  that  the  crack-tip  fields  are  captured,  and  large 
enough  to  avoid  the  relatively  large  errors  in  the  finite  element  approximation  of  the  fracture 
process  zone.  These  findings  are  verified  here,  and  similar  calculations  are  made  for  the 
element-failure  method. 

The  meshes  used  in  the  node-release  convergence  study  contain  370  nodes  and  680 
elements  (coarse),  1419  nodes  and  2720  elements  (medium),  and  5557  nodes  and  10,880 
elements  (fine).  For  the  element-failure  study,  they  contain  728  nodes  and  1361  elements 
(coarse),  2862  nodes  and  5536  elements  (medium),  and  1 1,162  nodes  and  21,952  elements  (fine). 
The  medium  meshes  are  those  of  Figure  6,  and  their  element  edge  lengths  are  approximately  one- 
half  those  of  the  coarse  meshes.  Likewise,  the  element  edge  lengths  of  the  fine  meshes  are  one- 
half  those  of  the  medium  meshes. 

The  top  of  Figure  10  is  a  plot  of  T*  versus  time  for  each  of  the  three  meshes  using  the 
node-release  method  with  a  zero  radius  of  T*.  The  bottom  of  Figure  10  is  the  same  plot  with  a 
radius  of  T*  of  0.5  cm.  Convergence  is  evident  when  the  radius  is  nonzero,  but  lacking  when  it 
is  zero.  These  results  verify  the  need  to  use  a  finite  value  of  the  radius  in  computations.  The 
limit  value  of  zero  in  the  definition  of  Equation  1  renders  nonconvergent  values  of  T*. 

The  same  comparison  can  be  made  for  the  element-failure  algorithm.  The  top  of  Figure 
1 1  is  a  plot  of  T*  versus  time  for  each  of  the  three  meshes,  using  the  element-failure  method  with 
a  zero  radius  of  T*.  The  bottom  of  Figure  11  is  the  same  plot  using  a  radius  of  0.5  cm.  As  was 
the  case  for  the  node-release  method,  T*  converges  for  the  nonzero  value  of  the  radius  of  T*,  but 
does  not  when  the  radius  is  zero.  Similar  behavior  is  obtained  from  the  node-release  and 
element-failure  methods,  indicating  that  the  cause  of  divergence  for  a  zero  radius  of  T*  is  not  a 
characteristic  of  the  crack-propagation  routine  (node-release  or  element-failure). 
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T*  (kN/m)  T*  (kN/m) 


Figure  10.  T  Versus  Time  with  Mesh  Refinement  for  Node  Release  and  Radius  of 
P  of  Zero  (top)  and  0.5  cm  (bottom)  iioa*_Tio3794  .u»] 
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T*  (kN/m)  T*  (kN/m) 


Figure  11. 


T*  Versus  Time  with  Mesh  Refinement  for  Element  Failure  and  Radius 
of  r*  of  Zero  (top)  and  0.5  cm  (bottom)  [iia/b_Tio3794  «f] 
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While  the  definition  of  T*  in  Equation  1  includes  the  limit  as  F*  goes  to  zero,  it  is 
apparent  by  the  demonstrated  lack  of  convergence  that  this  limit  cannot  be  achieved  by  simply 
using  a  zero  radius  of  T*  in  finite-element  modeling.  It  therefore  remains  to  choose  an  optimal 
nonzero  value  of  this  radius.  Since  the  energy  absorbed  by  the  plastic  deformation  of  elements  in 
the  exclusion  zone  is  erroneously  attributed  to  the  crack  tip  (for  the  formation  of  new  crack  faces) 
and  included  in  the  calculated  value  of  T*,  it  is  desirable  to  minimize  the  nonzero  radius  of  T*. 
Therefore,  a  recommended  value  of  this  radius  is  that  which  excludes  from  Cl  only  the  first  layer 
of  elements  around  the  crack  tip.  Additionally,  the  size  of  the  elements  around  the  crack  tip 
should  be  small  in  relation  to  the  size  of  the  plastic  zone. 

To  determine  the  sensitivities  of  the  T*  calculations  in  the  coupled  pressure-bar  specimen 
of  this  section,  the  medium  mesh  on  the  right  side  of  Figure  6  is  used  in  a  parametric  study  of  the 
radii  of  T  and  T*.  The  top  of  Figure  12  shows  plots  of  T*  versus  time  for  three  values  of  the 
radius  of  T*:  0.25,  0.50,  and  0.75  cm.  (The  radius  of  T  is  1.50  cm.)  The  value  of  0.25  cm  is 
chosen  to  represent  the  smallest  value  that  effectively  differs  from  zero,  inasmuch  as  it  excludes 
the  first  layer  of  elements  next  to  the  crack  tip.  The  value  of  0.50  cm  is  that  used  in  the  previous 
convergence  study.  The  value  of  0.75  cm  excludes  three  to  four  layers  of  elements  around  the 
crack  tip.  The  differences  in  these  plots  are  an  indication  of  the  size  of  the  complimentary  errors 
due  to  the  inability  of  the  mesh  to  resolve  the  crack  tip  fields  and  the  erroneous  contribution  of 
plastic  work  in  T*  to  the  crack  tip.  By  comparison,  the  bottom  of  Figure  12  shows  plots  of  T* 
versus  time  for  three  values  of  the  radius  of  T:  1.0,  1.5,  and  2.0  cm.  (The  radius  of  T*  is  0.50 
cm.)  The  value  of  1.0  cm  ranges  from  four  to  six  layers  of  elements  around  the  crack  tip.  Since 
there  is  no  reference  to  T  in  the  definition  of  T*,  valid  numerical  approximations  should  not  be 
sensitive  to  its  size.  For  this  example  with  substantial  crack-tip  plasticity,  the  calculated  values 
of  T*  appear  relatively  insensitive  to  changes  in  the  size  of  both  T  and  T*. 

2.3.3  Off-centered  Notch 

The  advantages  of  the  element-failure  method  lie  in  the  ease  with  which  it  models  crack 
propagation  in  general  directions.  In  this  section,  the  coupled  pressure-bar  specimen  of 
subsection  2.3.1  is  modified  so  that  the  notch  is  no  longer  centered  between  the  grips.  The  mesh, 
shown  at  the  top  of  Figure  13,  contains  the  same  number  of  nodes  and  elements  as  the  medium 
mesh  of  subsection  2.3.2.  The  loading  and  material  characterization  are  unchanged. 
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T*  (kN/m)  T*  (kN/m) 


Figure  12.  T*  Versus  Time  for  Three  Radii  of  V  (Top)  and  Three  Radii  of  r 

(Bottom)  [figl 2top/bottom.T 1 03794.tif) 
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Figure  13.  Off-centered  Notch  Meshes  at  0  |js  (top)  and  250  ps  (bottom,  failed 
elements  shaded)  [i2aA_Ti03794.tif] 

The  bottom  of  Figure  13  displays  the  mesh  at  250  ps,  with  the  failed  elements  shaded. 
Here  the  direction  of  propagation  was  determined  by  the  direction  of  maximum  T*.  These 
calculations  were  also  performed  with  the  direction  of  propagation  perpendicular  to  the 
maximum  normal  stress  in  the  region  of  the  crack  tip,  again  resulting  in  a  crack  path  that  veers 
toward  the  free  surface  on  the  right  side  of  the  specimen.  These  calculations  demonstrate  the 
ability  of  the  element-failure  method  (including  the  direction-determination  criterion)  to  model 
the  propagation  of  curved  cracks. 
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2.3.4  Interior  Crack 


The  preceding  calculations  of  coupled  pressure-bar  specimens  rely  on  a  notch  to  produce 
stress  concentrations  from  which  cracks  are  generated.  The  element-failure  method  also  provides 
a  convenient  method  for  modeling  interior  cracks  that  are  generated  from  a  material  flaw.  The 
flaw  can  be  simulated  by  failing  an  element  prior  to  loading.  Multiple  crack  tips  can  then  be 
allowed  to  propagate  from  the  flaw,  as  long  as  the  domain  Q,  used  to  calculate  T*  for  one  tip, 
does  not  capture  another  tip.  Figure  14  demonstrates  how  the  contours  can  be  drawn  for  the  two 
tips  of  an  interior  crack.  Note  that  the  domains  may  overlap. 


Q  of  lower  tip 

Figure  14.  Domains  of  Integration  for  Interior  Crack  [i 3_ti 03794 

A  sample  calculation  of  an  interior  crack  has  been  performed  on  a  rectangular  block  of 
4982  nodes  and  9760  elements.  The  block  measures  10  by  15  cm  and  the  material  model  is 
identical  to  that  of  subsection  2.3.1.  The  block  is  put  in  tension  by  applying  the  loading  of 
Equation  3  to  both  ends. 

Figure  15  shows  the  meshes  at  150  ps  (top)  and  200  ps  (bottom).  The  initial  flaw  is 
modeled  by  one  failed  element  adjacent  to  the  centerline.  Two  crack  tips  subsequently  propagate 
in  opposite  directions  (up  and  down)  perpendicular  to  the  tensile  stress  field,  as  determined  by 
the  direction  of  maximum  T*.  Due  to  symmetry,  an  exact  solution  would  indicate  identical 
values  of  T*  and  tip  velocity  for  both  crack  tips.  The  finite  element  solution  of  this 
demonstration  produces  the  values  of  T*  for  the  two  tips  plotted  in  Figure  16,  with  the  degree  of 
agreement  suggesting  the  ability  of  the  element-failure  method  to  model  interior  cracks. 
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Figure  15.  Tensile  Specimen  Meshes  at  150  ps  (top)  and  200  ps  (bottom) 

[14a/b_T1 03794.tif] 
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T*  (kN/rn) 


time  (s) 


Figure  16.  T*  Versus  Time  for  Top  and  Bottom  Tips  of  Interior  Crack  [i5_Ti03794.tif] 
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SECTION  3 
3D  ALGORITHM 


This  section  describes  an  extension  of  the  2D  algorithm  to  thin  structures  in  three 
dimensions.  Computations  are  then  presented  to  compare  the  accuracy  of  the  algorithm  to  a 
node-splitting  technique.  Also,  the  fracture  of  a  thin-walled  cylinder  loaded  in  torsion  is 
simulated  to  demonstrate  the  three-dimensional  capability. 

3.1  THE  3D  ELEMENT-FAILURE  ALGORITHM 

Figure  17  demonstrates  methods  of  modeling  cracks  with  finite  elements.  The  top  left  of 
the  figure  depicts  a  node-splitting  representation  of  a  crack  in  two  dimensions.  For  this  type  of 
representation,  the  newly  formed  crack  faces  are  modeled  by  contact  surfaces,  necessitating 
remeshing  around  the  crack  tip. 

The  top  right  of  Figure  17  depicts  a  crack  represented  by  the  element-failure  technique  in 
two  dimensions,  as  described  in  Section  2.  The  crack  tip  is  tracked  through  the  mesh,  and 
elements  are  failed  as  the  tip  passes  through  them.  Failed  elements  support  no  deviatoric  or 
tensile  volumetric  stresses,  modeling  frictionless  crack  faces  along  their  boundaries  with  unfailed 
elements.  It  can  be  seen  in  Figure  16  that  each  node  remains  attached  to  at  least  one  unfailed 
element  as  the  crack  tip  moves  through  the  mesh. 

The  extension  of  this  element-failure  algorithm  to  three  dimensions  is  depicted  at  the 
bottom  of  Figure  17.  The  extension  is  not  fully  general;  it  is  assumed  that  the  structure  contains 
one  small  thickness  relative  to  the  other  dimensions,  and  that  the  crack  front  spans  this  thickness. 
These  restrictions  logically  permit  the  computation  of  one  value  (averaged  over  the  length  of  the 
crack  front)  of  the  fracture  energy  per  unit  area  of  new  crack  face,  T*,  and  one  direction  of 
propagation.  All  points  along  the  crack  front  are  therefore  restricted  to  propagation  at  the  same 
rate  in  the  same  direction.  The  structural  thickness,  however,  need  not  be  uniform. 
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Figure  17.  Methods  of  Modeling  Cracks  with  Finite  Elements:  2D  Node  Splitting 

(top  left),  2D  Element  Failure  (top  right)  and  3D  Element  Failure  (bottom, 
unfailed  elements  not  shown) 
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To  most  accurately  model  a  crack  front  of  uniform  speed  and  direction  of  propagation, 
the  elements  must  be  aligned  in  columns  through  the  thickness  of  the  structure,  as  shown  at  the 
bottom  of  Figure  17.  A  direct  benefit  of  this  element  arrangement  is  that  the  element-failure 
algorithm  introduced  for  planar  models  can  be  applied  to  the  surface  of  the  3D  structure  in  a 
straightforward  manner.  The  extension  to  three  dimensions  is  then  reduced  to  including  the 
effects  of  the  elements  that  are  aligned  with  the  top-surface  elements  in  the  routines  applied  to 
the  top  surface.  Such  an  implementation  will  be  described  in  the  following  two  subsections. 


Under  the  current  restrictions  (a  thin  3D  structure  with  a  crack  front  of  uniform  speed  and 
direction  of  propagation  spanning  the  thin  dimension),  it  would  also  be  appropriate  to  model  the 
3D  structure  with  shell  elements  and  apply  the  2D  element-failure  method  directly  to  the  shell 
elements.  The  approach  of  using  solid  elements  aligned  in  columns  through  the  thin  dimension 
of  the  structure  is  used  here  due  to  the  ease  with  which  the  solid  elements  can  be  attached  to 
other  parts  of  the  structure  which  do  not  possess  a  thin  dimension  and  must  be  modeled  by  solid 
elements.  An  example  of  such  a  geometry  is  that  of  the  steel  case  of  a  hard  target  penetrator. 

3.1 .1  Fracture  Model 

The  bottom  of  Figure  17  demonstrates  a  trail  of  failed  elements  traversed  by  a  crack  front, 
and  how  the  failed  elements  model  a  crack.  To  simulate  the  flaw  required  to  initiate  fracture,  the 
algorithm  fails  a  column  of  elements  specified  by  the  user.  If  this  flaw  lies  on  an  edge  of  the 
structure,  one  value  of  T*  is  calculated  and  a  corresponding  crack  front  may  propagate  through 
the  structure,  simulating  an  edge  crack.  If  the  flaw  is  not  on  an  edge,  two  values  of  T  are 
calculated  and  two  corresponding  crack  fronts  may  propagate,  simulating  an  internal  crack. 


The  “top”  surface  of  the  structure  denotes  the  surface  on  which  the  planar  element-failure 
algorithm  is  applied.  The  intersection  of  the  crack  front  and  the  top  surface  is  used  to  define  the 
crack  path,  as  in  Figure  17.  The  path  is  restricted  to  travel  in  a  straight  line  across  each  element 
face.  On  the  timestep  during  which  it  enters  a  new  element  face,  a  routine  to  determine  the  new 
direction  of  propagation  is  invoked.  In  this  section,  the  direction  of  propagation  is  taken  as  the 
direction  normal  to  the  maximum  principal  stress  (mode  I  fracture)  for  a  state  of  stress  that  is 
obtained  by  averaging  the  stresses  along  the  crack  front.  While  the  direction  of  propagation  is 
only  updated  when  the  crack  front  enters  a  new  column  of  elements,  the  fracture  energy  is 
calculated  at  each  timestep,  and  —  combined  with  the  fracture  criterion  it  updates  the  speed  of 
propagation. 
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Although  the  crack  path  is  defined  on  the  top  surface  of  the  structure,  it  represents  the 
history  of  the  entire  crack  front.  Due  to  the  alignment  of  elements  in  columns  through  the 
thickness,  when  the  crack  path  enters  a  new  element  face,  the  crack  front  simultaneously  enters  a 
new  column  of  elements.  An  important  consequence  is  that  all  elements  through  the  thickness 
can  be  failed  simultaneously,  maintaining  a  good  approximation  to  a  crack  front  of  uniform 
speed  and  direction  of  propagation. 

The  failed  elements  are  not  permitted  to  support  deviatoric  or  tensile  volumetric  stresses, 
and  they  can  be  realized  by  simply  modifying  their  constitutive  relations  so  that  they  can  sustain 
only  positive  pressure,  typically  occurring  under  negative  volumetric  strains.  When  an  element  is 
failed,  its  stresses  are  set  to  zero  and  their  effect  is  replaced  by  externally  applied  nodal  forces 
that  are  equivalent  to  the  internal  nodal  forces  contributed  by  the  element  stresses  before  zeroing. 
These  external  nodal  forces  are  reduced  to  zero  linearly  as  the  crack  front  traverses  the  element. 
They  reduce  the  release  waves  associated  with  the  sudden  failure  of  elements  [7,  8,  9]  and  extract 
energy  from  the  structure.  A  surface  between  a  failed  and  an  unfailed  element  acts  as  a  traction- 
free  crack  face  when  the  volumetric  strain  of  the  failed  element  is  positive  (opened  crack),  and  a 
pressure-transmitting  crack  face  when  the  volumetric  strain  of  the  failed  element  is  negative 
(closed  crack). 

3.1.2  The  T*  Fracture  Parameter 

As  a  measure  of  the  energy  release  associated  with  a  unit  area  of  crack  surface  extension 
in  dynamic  fracture,  T*  [17, 18]  is  expressed  as, 

T*  =  l  [(w  +  T)n„  -  <v,u  Jdr  *  (4) 

which  is  identical  to  Equation  1  for  2D  geometry,  except  that  T*  is  now  a  surface  contour  around 
the  crack  front.  W  and  T  represent  the  internal  energy  density  and  kinetic  energy  density, 
respectively;  ay,  the  ij-th  component  of  the  stress  tensor;  nj,  the  j-th  component  of  the  unit 
normal  vector  to  T*;  np,  the  component  of  the  unit  normal  vector  in  the  direction  of  crack 
propagation;  Uj,  the  i-th  component  of  the  displacement  vector;  uyp,  the  derivative  of  the  i-th 
component  of  the  displacement  vector  in  the  direction  of  crack  propagation. 
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The  largest  errors  in  a  numerical  solution  occur  in  the  region  of  the  crack  front  where 
field  quantities  vary  most  rapidly.  As  a  result,  it  is  again  desirable  to  equate  the  limiting  integral 
of  Equation  4  to  volume  and  surface  integrals  by  the  divergence  theorem  [7], 

T*  =  Jr[(W  +  T)np  -ajjnjUj  JdT  -  ^[(w  +  T^p-ptijUjp -c^U; j),p]dfl  (5) 

where  T  represents  a  far-field  crack-front  surface  contour  and  Q  is  the  volume  between  T  and  T*. 
Here,  nj  is  the  j-th  component  of  the  unit  normal  vector  to  T.  (A  subscript  following  a  comma 
represents  the  direction  of  a  spatial  derivative.) 

The  top  of  Figure  18  shows  an  example  of  the  domains  of  T,  T*  and  Q  for  the  calculation 
of  T*  in  a  planar  model,  and  the  bottom  of  Figure  17  shows  them  (without  elements)  about  a 
crack  front  spanning  the  thin  dimension  of  a  structure.  Although  T*  is  defined  over  T*  in  the 
limit  as  it  is  shrunk  about  the  crack  front  (Equation  4),  only  the  domains  of  T  and  Q  are  used  in 
the  computations  (Equation  5). 

To  determine  which  elements  comprise  Q  and  which  element  faces  comprise  T,  the 
algorithm  developed  for  the  planar  models  can  be  applied  to  the  element  faces  of  the  top  surface. 
In  this  way,  the  distance  R  between  the  centroid  of  each  top-surface  element  face  and  the  crack 
front  at  the  top  surface  is  compared  to  nominal  values  of  the  radius  of  T*,  Rr*>  and  the  radius  of 
T,  Rr-  If  R  lies  between  the  two,  Rr*  ^  R  ^  Rr>  then  the  top-surface  element  and  all  elements 
aligned  with  it  are  included  in  Q.  Figure  19  depicts  this  comparison  for  one  top-surface  element 
and  the  elements  aligned  with  it. 

Once  the  comparison  has  been  made  for  all  top-surface  elements  and  the  domain  Q  has 
therefore  been  established,  the  element  faces  of  Q  are  then  categorized  as  either  members  of  T* 
or  T.  Those  that  define  the  crack-front  surface  contour  belong  to  T*;  all  others  belong  to  T.  It 
has  been  recommended  [7,  8]  that  all  elements  that  have  been  inside  F*  be  excluded  from  Q  in 
subsequent  evaluations  of  T*  due  to  the  lack  of  accuracy  of  the  numerical  solution  in  the  crack- 
front  region  and  the  continuation  of  those  errors  when  path-dependent  materials  (plasticity)  are 
unloaded.  The  result  is  a  zone  of  elements  that  are  not  included  in  Q  forming  a  wake  behind  the 
crack  front.  The  boundary  of  this  exclusion  zone  defines  T*. 
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Figure  19.  Distance  from  Crack  Front  for  Candidate  Elements  of  Q 

Across  each  tetrahedral  element,  the  displacement  is  interpolated  linearly,  resulting  in 
constant  displacement  derivatives,  stresses,  densities,  and  internal  energy  densities.  Numerical 
integration  of  the  integral  over  T  in  Equation  5  is  therefore  achieved  by  one-point  integration  on 
each  of  the  boundary  surfaces. 

The  integral  over  Q  in  Equation  5  is  integrated  by  employing  one-point  quadrature  over 
each  tetrahedral  element  in  Q.  For  constant-strain  tetrahedra,  the  second-order  displacement 
derivatives  must  be  smoothed  before  application  of  the  quadrature  scheme  (because  they  are  zero 
in  the  interiors  of  the  elements  and  infinite  at  the  element  boundaries).  Smoothing  is  achieved  by 
averaging  values  of  the  first  displacement  derivatives  from  the  adjoining  elements  at  each  node, 
and  then  interpolating  those  nodal  values  linearly  over  each  element  to  render  constant  second 
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displacement  derivatives  in  each  element.  A  similar  approach  is  used  for  the  energy-density 
derivatives. 

It  has  been  demonstrated  [8]  that  a  lack  of  convergence  with  mesh  refinement  occurs  for 
solutions  of  ductile  fracture  if  Rr*  is  taken  as  the  limit  value  of  zero.  For  this  reason,  the 
computations  of  the  following  section  were  performed  with  a  value  of  Rr*  just  large  enough  to 
exclude  one  or  two  columns  of  elements  on  all  sides  of  the  crack  front  from  £2. 

3.2  NUMERICAL  RESULTS 

As  noted  previously,  the  3D  algorithm  has  been  incorporated  into  the  1999  version  of  the 
EPIC  code  [2,  19].  In  the  example  computations  that  follow,  the  direction  of  propagation  is 
taken  normal  to  the  maximum  principal  stress. 

3.2.1  Comparison  of  3D  Element  Failure  and  2D  Node  Splitting 

The  primary  advantage  of  the  element-failure  method  is  the  simplicity  of  implementation. 
This  subsection  is  intended  to  demonstrate  its  accuracy.  It  has  been  shown  in  Section  2,  by  the 
comparison  of  coupled  pressure-bar  simulations,  that  the  accuracy  of  the  element-failure  method 
presented  for  planar  problems  is  similar  to  that  of  the  node-splitting  method.  Here,  the  same 
coupled  pressure  bar  will  be  modeled  in  three  dimensions  and  results  will  be  compared  to  plane- 
strain  node-splitting  results. 

The  coupled  pressure  bar  experiment  is  depicted  in  Figure  20.  The  geometry  is  identical 
to  that  used  for  the  2D  computations  in  Section  2,  except  that  the  3D  geometry  requires  a 
specified  thickness.  The  material  model  and  loading  are  identical  to  those  used  for  the  2D 
problem  in  Section  2. 

The  mesh  of  the  node-splitting  algorithm  is  shown  on  the  left  side  of  Figure  21,  and  it 
contains  1419  nodes  and  2720  elements.  Symmetry  allows  the  modeling  of  half  the  specimen, 
and  node  “splitting”  is  achieved  by  releasing  nodes  along  the  line  of  symmetry,  replacing  each 
restraint  by  an  external  nodal  force  that  is  equivalent  to  the  internal  nodal  force  due  to  element 
stresses,  and  reducing  the  external  forces  to  zero  as  the  crack  tip  traverses  the  element  edge  (as 
described  in  Section  2). 
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Figure  20.  Coupled  3D  Pressure-bar  Specimen  and  Fracture  Response  Function 


Figure  21.  Meshes  for  Comparison  of  2D  Node-splitting  (left,  2720  elements)  and 
3D  Element-failure  (right,  66,432  elements)  Methods 


For  the  element-failure  algorithm  in  three  dimensions,  no  symmetry  assumptions  are 
made  and  a  mesh  of  17,076  nodes  and  66,432  elements  is  used,  as  shown  on  the  right  side  of 
Figure  21.  The  mesh  contains  5536  triangular  element  faces  on  the  top  surface,  as  compared  to 
the  5440  triangular  elements  represented  by  the  plane-strain  node-splitting  mesh  (left  side  of 
Figure  21). 

The  nodes  are  restrained  from  motion  normal  to  the  top  surface  to  simulate  plane  strain 
conditions.  For  both  algorithms,  Rp  =  1.5  cm  and  Rr*  =  0.5  cm. 
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Figure  22  contains  plots  of  T*  and  the  z-coordinate  of  the  crack  front  as  functions  of  time 
for  both  the  node-splitting  and  element-failure  methods.  Similar  accuracy  of  the  two  methods  is 
reflected  in  the  close  agreement  of  the  solutions. 
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Figure  22.  Comparison  of  Results  for  2D  Node-splitting  and  3D  Element-failure 
Methods 
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Figure  23  is  a  plot  of  the  element-failure  mesh  at  250  ps,  with  the  failed  elements 
darkened.  The  crack  can  be  seen  to  propagate  up  the  line  of  symmetry  with  little  variation.  , 


Figure  23.  Element-failure  Mesh  at  250  |is  (failed  elements  darkened) 

3.2.2  Fracture  of  Thin-walled  Cylinder  under  Torsion 

The  computation  in  this  subsection  is  intended  to  demonstrate  the  three-dimensional 
capability  of  the  algorithm.  The  problem  consists  of  a  thin-walled  cylinder  under  torsional 
loading.  One  end  of  the  cylinder  is  restrained  against  motion  in  the  plane  normal  to  the  cylinder 
axis,  while  the  other  end  is  subjected  to  a  rotational  velocity,  as  shown  in  Figure  24.  Both  ends 
are  free  to  deform  in  the  axial  direction. 

The  elastic/perfectly  plastic  material  model  of  the  previous  subsection  is  also  used  for 
these  computations.  The  mesh  is  composed  of  20,376  nodes  and  69,120  tetrahedral  elements,  as 
shown  on  the  left  side  of  Figure  25.  An  initial  flaw  in  the  material  is  simulated  by  failing  a 
column  of  elements  near  the  midlength  of  the  cylinder.  This  column  is  darkened  in  the  plot  on 
the  left  side  of  Figure  25.  In  the  computations  of  T*,  Rr  =  1.5  cm  and  Rr*  =  0.5  cm. 
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Figure  25.  Thin-walled  Cylinder  Mesh  at  0  ps  (left)  325  ps  (center),  and  400  ps 

(right).  Failed  elements  are  darkened,  including  the  initial  flaw  at  left. 
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After  325  jus  of  loading,  two  crack  tips  have  begun  to  propagate  in  opposite  directions 
from  the  flaw.  The  plot  at  the  center  of  Figure  25  shows  the  mesh  at  this  time,  with  failed 
elements  darkened.  After  400  ps,  the  crack  tips  have  propagated  farther,  as  shown  on  the  right 
side  of  Figure  25,  and  the  crack  is  opening  in  mode  I.  Only  those  elements  through  which  the 
crack  front  has  passed  have  been  failed.  The  increased  width  of  the  crack  at  400  ps  is 
represented  by  the  failed  elements  that  are  allowed  to  freely  deform  under  the  zero  stress  state 
associated  with  expanded  failed  elements. 
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SECTION  4 

PARAMETRIC  COMPUTATIONAL  STUDY 


The  final  objective  for  this  work  is  to  have  the  capability  to  perform  3D  crack 
propagation  computations  for  penetrators  impacting  a  variety  of  hard  concrete  targets.  The  series 
of  parametric  computations  that  follows  demonstrates  this  capability. 

A  description  of  the  configuration  used  for  the  parametric  computations  is  shown  m 
Figure  26.  The  initial  flaw  positions  occur  at  points  A,  B,  C,  and  D.  The  geometry  and  impact 
conditions  of  the  generic  hard  target  penetrator  are  similar  to  that  presented  previously  by 
Johnson,  et.  al.  [20].  The  penetrator  has  a  length  of  L  =  203  cm  and  a  diameter  of  D  =  25.4  cm, 
for  L/D  =  8.  The  length  of  the  nose  is  63.5  cm  and  the  wall  thickness  is  2.54  cm.  The  case 
material  is  4340  steel  with  a  hardness  of  Rc49.  A  characterization  of  this  material  for  the 
Johnson-Cook  strength  model  is  provided  in  [21],  and  it  is  represented  by  material  number  85  in 
the  EPIC  material  library.  The  inert  explosive  is  Filler  E,  represented  by  material  number  60  in 
the  library.  The  total  weight  is  366  kg. 

The  penetration  is  represented  by  65,736  tetrahedral  elements.  The  inert  explosive  and 
the  penetrator  case  utilize  the  elements  in  a  symmetric  arrangement  (24  tets  per  brick),  and  the 
steel  nose  uses  a  nonsymmetric  arrangement  (6  tets  per  brick).  Symmetric  bricks  were  used  for 
the  case  because  it  was  determined  that  the  crack  propagation  computations  were  sometimes 
more  accurate  and  robust  with  this  arrangement.  The  newly  developed  uniform  ring  option  for 
rod  shapes  (in  the  EPIC  preprocessor)  was  used  to  meet  the  through  thickness  assumptions  in  the 
algorithm.  The  nonsymmetric  brick  arrangement  was  used  in  the  nose  because  there  is  no  current 
capability  to  generate  symmetric  bricks  with  a  uniform  ring  arrangement  for  nose  shapes.  Also, 
the  entire  penetrator  (without  a  plane  of  symmetry)  is  included  in  the  model  because  the  crack 
propagation  algorithm  cannot  accurately  account  for  planes  of  symmetry.  Furthermore,  after  a 
crack  begins  to  form  there  is  no  longer  a  plane  of  symmetry. 

The  crack  propagation  model  (crack  propagation  velocity  versus  T*)  is  discussed  and 
presented  in  Part  1  of  this  report  [1].  The  relationship  between  T*  and  crack  velocity  is  also 
included  in  Figure  26. 
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The  concrete  target  has  a  compressive  strength  of  fc'  =  34.5  MPa  (5000  lb/in2),  and  a 

density  of  2280  kg/m3.  It  is  assumed  to  be  infinitely  thick  and  infinitely  wide.  The  resistance  of 
the  target  is  represented  by  the  PENCRV3D  model  that  is  linked  to  the  finite  element 
representation  of  the  penetrator  [2].  This  simplified  representation  of  the  target  allows  the 
computations  to  be  run  much  faster  (hours  instead  of  days  of  CPU  time),  such  that  parametric 
computations  can  be  performed. 

Figure  27  shows  the  dynamic  response  of  the  baseline  computation,  in  which  cracks  are 
excluded.  At  2.0  ms  the  penetrator  is  loaded  such  that  bending  introduces  compression  in  the 
forward  wall  of  the  penetrator  and  tension  in  the  aft  wall.  At  5.0  ms  the  penetrator  straightens  to 
approximate  the  original  configuration,  and  at  10.0  ms  it  has  almost  stopped. 

Figure  28  shows  the  computed  response  for  the  initial  flaw  at  point  A.  At  a  time  of 
approximately  1.5  ms,  T*  begins  to  rise  for  both  fronts  and  the  crack  propagates.  At  about 
2.3  ms  T*  for  both  fronts  goes  negative  and  the  crack  no  longer  propagates.  At  2.8  ms  T*  for  the 
aft  front  goes  positive  and  the  computation  stops  because  the  prescribed  array  size  for  the  number 
of  cracked  elements  is  exceeded. 

Figure  29  shows  the  computed  response  for  a  flaw  at  point  B.  Here  the  crack  propagates 
under  the  positive  T*  that  occurs  between  1.5  and  2.0  ms.  When  the  aft  crack  tip  reaches  the 
base  plate  of  the  penetrator  it  can  no  longer  propagate.  The  computation  continues,  however,  out 
to  the  full  10.0  ms. 

Figure  30  shows  the  response  for  a  flaw  at  point  C.  Here  a  large  value  of  T*  is  computed 
and  the  crack  propagates  rapidly.  For  this  condition  the  computation  is  discontinued  at  2.1  ms. 
A  probable  explanation  is  that  the  reversal  in  loading  tends  to  change  the  direction  of  the  crack 
such  that  it  tries  to  propagate  back  into  the  previously  cracked  elements.  The  algorithm  in  the 
code  cannot  currently  handle  this  condition,  and  the  computation  is  discontinued  when  this 
occurs. 


Figure  31  shows  the  response  for  a  flaw  at  point  D.  The  location  of  this  flaw  is 
essentially  under  compression  and  does  not  propagate  significantly  because  T*  is  generally 
negative.  Again,  it  is  probable  that  the  reversal  of  loading  is  what  caused  this  computation  to 
stop  at  4.7  ms.  It  is  interesting  to  compare  and  contrast  the  propagation  at  point  C  (extensive 
propagation  in  a  tensile  field)  to  that  at  point  D  (little  propagation  in  a  compressive  field). 
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It  is  difficult  to  assess  the  accuracy  of  these  computations  because  there  are  no  known  test 
data  and/or  analytic/computational  solutions  for  curved  surfaces  under  reversing  loading 
conditions.  The  computed  values  of  T*  are  generally  very  high  and  produce  more  crack 
propagation  than  would  be  expected. 

Several  questions  that  need  to  be  answered  are  as  follows: 

•  What  is  the  effect  of  curvature. 

•  What  is  the  effect  of  zoning. 

.  What  is  the  effect  of  general  loading  that  may  reverse  (and  is  not  always  increasing). 

•  What  is  the  effect  of  through  thickness  variations  in  stress  and  strain. 

•  What  is  the  effect  of  the  ratio  of  the  case  radius  to  the  case  wall  thickness. 

•  What  is  the  effect  of  the  ratio  of  the  case  radius  to  the  T*  outer  radius,  T. 

Future  work  should  address  these  issues  such  that  the  accuracy  and/or  robustness  of  the 
algorithm  can  be  better  evaluated  and/or  improved  for  hard  target  penetrator  applications. 
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Figure  26.  Description  of  the  Parametric  Computations 
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forward  tip 


Cross-section  External  View 


Figure  31.  Computation  for  Initial  Flaw  at  Point  D 
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SECTION  5 

SUMMARY  AND  CONCLUSIONS 

This  report  has  provided  2D  and  3D  crack  propagation  algorithms  for  general  directions, 
as  well  as  parameteric  crack  propagation  computations  for  a  penetrator  impacting  a  hard  concrete 

target. 

The  element-failure  method  of  modeling  crack  propagation  assimilates  the  fracture 
mechanics  of  node-splitting  methods  with  the  simplicity  of  micromechanically  based  methods. 
The  T*  fracture  parameter,  a  measure  of  the  energy  release  associated  with  the  creation  of  a  unit 
area  of  fracture  surface  under  dynamic  loading  represents  a  fracture  mechanics  approach  to  crack 
propagation.  The  technique  of  failing  elements  eliminates  the  need  to  remesh  around  the  crack 
tip  and  redefine  contact  surfaces. 

Numerical  examples  show  similar  accuracy  of  the  node-splitting  (node-release,  along  a 
line  of  symmetry)  method  and  the  element-failure  method  for  meshes  of  similar  refinement. 
Convergence  studies  verify  the  need  for  a  nonzero  size  of  the  contour  T*,  used  in  the  evaluation 
of  T*,  for  both  methods.  Convergence  with  mesh  refinement  is  lost  if  the  size  is  zero.  The 
direction  of  propagation  of  the  element-failure  method  is  chosen  as  the  direction  normal  to  the 
maximum  stress,  or  the  direction  of  the  maximum  value  of  T*.  The  ability  of  the  element-failure 
method  to  model  interior  cracks  is  demonstrated  assuming  the  existence  of  an  initial  element¬ 
sized  material  flaw. 

The  most  important  characteristic  of  the  element-failure  method  is  that  it  allows  for  three 
dimensional  computations  to  be  performed  without  rezoning.  The  3D  algorithm  has  been  shown 
to  provide  good  agreement  with  2D  computations,  and  to  provide  the  capability  to  propagate 
cracks  in  general  directions.  Although  the  3D  algorithm  has  provided  accurate  and  robust 
behavior  for  flat  plates,  the  application  to  curved  surfaces  and  changing  loading  directions  has 
raised  some  questions  that  have  not  yet  been  fully  resolved. 

Finally,  the  3D  algorithm  has  been  successfully  applied  tp  the  application  of  a  penetrator 
impacting  a  hard  concrete  target. 
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